#!/usr/bin/perl

use strict;
use PGPLOT;

# this script plots JR's cm model, RJS's K band model and JBS's mm model for 1934-638

# set up the array of frequencies
my $min_freq=1000; # in MHz
my $max_freq=100000; # in MHz

# do it in log space?
my $do_log=1;
my $write_file=0;

if ($do_log==1){
    $min_freq=log($min_freq)/log(10.0);
    $max_freq=log($max_freq)/log(10.0);
}

# how many points to make along the way
my $npoints=1000;
my @frequencies;
my @fluxes_JR;
my @fluxes_RJS;
my @fluxes_JBS;
my @fluxes_combined;
my @diff_JR;
my @diff_RJS;
my @diff_JBS;
my ($minimum_flux,$maximum_flux);
my ($minfraction,$maxfraction);
if ($write_file==1){
    open(OUT,">model_fluxes.txt");
}
for (my $i=0;$i<$npoints;$i++){
    $frequencies[$i]=$min_freq+$i*($max_freq-$min_freq)/$npoints;
    if ($do_log==1){
	$fluxes_JR[$i]=-30.7667+26.4908*$frequencies[$i]-
	    7.0977*($frequencies[$i]**2)+0.605334*($frequencies[$i]**3);
	$fluxes_RJS[$i]=-202.6259+149.7321*$frequencies[$i]-
	    36.4943*($frequencies[$i]**2)+2.9372*($frequencies[$i]**3);
	$fluxes_JBS[$i]=5.65274-1.31408*$frequencies[$i];
#	$fluxes_combined[$i]=-18.283+15.949*$frequencies[$i]-
#	    4.14656*($frequencies[$i]**2)+0.331408*($frequencies[$i]**3);
	$fluxes_combined[$i]=26.7035-15.6502*$frequencies[$i]+
	    3.23543*($frequencies[$i]**2)-0.242125*($frequencies[$i]**3);
#	$fluxes_JBS[$i]=$fluxes_combined[$i];
#	$fluxes_combined[$i]=-61.6339+61.5816*$frequencies[$i]-
#	    22.0199*($frequencies[$i]**2)+3.41902*($frequencies[$i]**3)-
#	    0.198525*($frequencies[$i]**4);
	$diff_JR[$i]=abs((10**($fluxes_JR[$i])-10**($fluxes_combined[$i]))/
			 10**($fluxes_JR[$i]));
	$diff_RJS[$i]=abs((10**($fluxes_RJS[$i])-10**($fluxes_combined[$i]))/
			  10**($fluxes_RJS[$i]));
	$diff_JBS[$i]=abs((10**($fluxes_JBS[$i])-10**($fluxes_combined[$i]))/
			  10**($fluxes_JBS[$i]));
    } else {
	$fluxes_JR[$i]=10**(-30.7667+26.4908*log($frequencies[$i])/log(10.0)-
			    7.0977*((log($frequencies[$i])/log(10.0))**2)+
			    0.605334*((log($frequencies[$i])/log(10.0))**3));
	$fluxes_RJS[$i]=10**(-202.6259+149.7321*log($frequencies[$i])/log(10.0)-
			     36.4943*((log($frequencies[$i])/log(10.0))**2)+
			     2.9372*((log($frequencies[$i])/log(10.0))**3));
	$fluxes_JBS[$i]=10**(5.65274-1.31408*log($frequencies[$i])/log(10.0));
# 	$fluxes_combined[$i]=10**(-18.283+15.949*log($frequencies[$i])/log(10.0)-
#				  4.14656*((log($frequencies[$i])/log(10.0))**2)+
#				  0.331408*((log($frequencies[$i])/log(10.0))**3));
 	$fluxes_combined[$i]=10**(26.7035-15.6502*log($frequencies[$i])/log(10.0)+
				  3.23543*((log($frequencies[$i])/log(10.0))**2)-
				  0.242125*((log($frequencies[$i])/log(10.0))**3));
#	$fluxes_combined[$i]=10**(-61.6339+61.5816*log($frequencies[$i])/log(10.0)-
#				  22.0199*((log($frequencies[$i])/log(10.0))**2)+
#				  3.41902*((log($frequencies[$i])/log(10.0))**3)-
#				  0.198525*((log($frequencies[$i])/log(10.0))**4));
	$diff_JR[$i]=abs(($fluxes_JR[$i]-$fluxes_combined[$i])/$fluxes_JR[$i]);
	$diff_RJS[$i]=abs(($fluxes_RJS[$i]-$fluxes_combined[$i])/$fluxes_RJS[$i]);
	$diff_JBS[$i]=abs(($fluxes_JBS[$i]-$fluxes_combined[$i])/$fluxes_JBS[$i]);
    }
    if ($write_file==1){
	if ((($do_log==0)&&($frequencies[$i]<10000))||
	    (($do_log==1)&&($frequencies[$i]<(log(10000.0)/log(10.0))))){
#	print OUT $frequencies[$i]." ".$fluxes_JR[$i]."\n";
	} elsif ((($do_log==0)&&($frequencies[$i]<25000))||
		 (($do_log==1)&&($frequencies[$i]<(log(25000)/log(10.0))))){
	    print OUT $frequencies[$i]." ".$fluxes_RJS[$i]."\n";
	} else {
	    print OUT $frequencies[$i]." ".$fluxes_JBS[$i]."\n";
	}
    }
    if ($i==0){
	$minimum_flux=$maximum_flux=$fluxes_JR[$i];
	$minfraction=$maxfraction=$diff_JR[$i];
    } 
    if ($fluxes_JR[$i]<$minimum_flux){
	$minimum_flux=$fluxes_JR[$i];
    }
    if ($fluxes_RJS[$i]<$minimum_flux){
	$minimum_flux=$fluxes_RJS[$i];
    }
    if ($fluxes_JBS[$i]<$minimum_flux){
	$minimum_flux=$fluxes_JBS[$i];
    }
    if ($fluxes_combined[$i]<$minimum_flux){
	$minimum_flux=$fluxes_combined[$i];
    }
    if ($fluxes_JR[$i]>$maximum_flux){
	$maximum_flux=$fluxes_JR[$i];
    }
    if ($fluxes_RJS[$i]>$maximum_flux){
	$maximum_flux=$fluxes_RJS[$i];
    }
    if ($fluxes_JBS[$i]>$maximum_flux){
	$maximum_flux=$fluxes_JBS[$i];
    }
    if ($fluxes_combined[$i]>$maximum_flux){
	$maximum_flux=$fluxes_combined[$i];
    }
    if ($diff_JR[$i]<$minfraction){
	$minfraction=$diff_JR[$i];
    }
    if ($diff_RJS[$i]<$minfraction){
	$minfraction=$diff_RJS[$i];
    }
    if ($diff_JBS[$i]<$minfraction){
	$minfraction=$diff_JBS[$i];
    }
    if ($diff_JR[$i]>$maxfraction){
	$maxfraction=$diff_JR[$i];
    }
    if ($diff_RJS[$i]>$maxfraction){
	$maxfraction=$diff_RJS[$i];
    }
    if ($diff_JBS[$i]>$maxfraction){
	$maxfraction=$diff_JBS[$i];
    }

}
$minfraction=-0.01;
$maxfraction=0.1;
if ($write_file==1){
    close(OUT);
}

# plot the models
#pgopen("44/xs");
#pgopen("1934_models_4th_5-25.ps/cps");
#pgopen("1934_models_3rd_1-100.ps/cps");
pgopen("1934_models_C007_new.ps/cps");
my ($vp_xmin,$vp_xmax,$vp_ymin,$vp_ymax);
pgqvp(0,$vp_xmin,$vp_xmax,$vp_ymin,$vp_ymax);
my ($cvp_ymin,$cvp_ymax);
$cvp_ymax=$vp_ymax;
$cvp_ymin=$vp_ymin+($vp_ymax-$vp_ymin)/5.0;
pgsvp($vp_xmin,$vp_xmax,$cvp_ymin,$cvp_ymax);
pgswin($min_freq,$max_freq,$minimum_flux,$maximum_flux);
pgsci(1);
if ($do_log==1){
    pgbox("BCLTS",0.0,0,"BCLNTS",0.0,0);
} else {
    pgbox("BCTS",0.0,0,"BCNTS",0.0,0);
}
pglab("","Flux (Jy)","1934-638 Flux Model Comparison");
pgsci(2);
pgsls(2);
pgline($npoints,\@frequencies,\@fluxes_JR);
pgsci(3);
pgline($npoints,\@frequencies,\@fluxes_RJS);
pgsci(5);
pgline($npoints,\@frequencies,\@fluxes_JBS);
pgsci(4);
pgsls(1);
pgline($npoints,\@frequencies,\@fluxes_combined);

my ($pglen_x,$pglen_y);
pglen(4,"combined model",$pglen_x,$pglen_y);
$pglen_y/=5.0;
print "length x = $pglen_x, y = $pglen_y\n";
my ($textloc_x,$textloc_y,$textloc_height);
$textloc_x=$max_freq-1.1*$pglen_x;
$textloc_height=1.1*$pglen_y;
$textloc_y=$maximum_flux-$textloc_height;
pgsci(2);
pgtext($textloc_x,$textloc_y,"JR model");
pgsci(3);
$textloc_y-=$textloc_height;
pgtext($textloc_x,$textloc_y,"RJS model");
#pgsci(4);
#$textloc_y-=$textloc_height;
#pgtext($textloc_x,$textloc_y,"JBS model");
pgsci(4);
$textloc_y-=$textloc_height;
#pgtext($textloc_x,$textloc_y,"combined model");
pgtext($textloc_x,$textloc_y,"JBS model");

$cvp_ymax=$cvp_ymin;
$cvp_ymin=$vp_ymin;
pgsvp($vp_xmin,$vp_xmax,$cvp_ymin,$cvp_ymax);
pgswin($min_freq,$max_freq,$minfraction,$maxfraction);
pgsci(1);
pgsch(0.8);
if ($do_log==1){
    pgbox("BCLNTS",0.0,0,"BCNTS",0.0,0);
} else {
    pgbox("BCNTS",0.0,0,"BCNTS",0.0,0);
}
pgsch(1.0);
pglab("Frequency (MHz)","Fractional Difference","");
pgsci(2);
pgline($npoints,\@frequencies,\@diff_JR);
pgsci(3);
pgline($npoints,\@frequencies,\@diff_RJS);
#pgsci(4);
pgsci(4);
pgline($npoints,\@frequencies,\@diff_JBS);

pgclos();
